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ABSTRACT 

The  problem  of  obtaining  a  control  policy  which  is  optimum  in  some 
respect  for  a  specified  system  is  one  of  the  more  pressing  problems  of 
automatic  control  theory  today.   L.  S.  Pontryagin  stated  his  Maximum 
Principle  in  1956  and  L.  I.  Rozonoer  discussed  and  extended  the  work  in 
a  series  of  articles  published  in  1959. 

The  Maximum  Principle,  as  it  applies  to  the  "free  right  end"  problem 
for  a  non-linear,  one-dimensional  system,  is  utilized  in  conjunction  with 
the  CDC  1604  digital  computer  to  obtain  optimal  control  policies  for 
several  different  cost  functions. 

The  authors  wish  to  express  their  appreciation  for  the  encourage- 
ment given  them  by  Dr.  Harold  A.  Titus  of  the  U.  S.  Naval  Postgraduate 
School  and  to  express  their  gratitude  to  Miss  Mary  E.  Haynes  of  the 
Computer  Facility,  U.  S.  Naval  Postgraduate  School  whose  patient  program- 
ming assistance  materially  aided  in  the  successful  completion  of  this 
thesis. 
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1.   Introduction. 

One  of  the  most  important  problems  of  automatic  control  theory  is 
the  problem  of  creating  systems  which  are  optimal  in  some  prescribed 
sense.   Since  the  advent  of  satellites  and  space  travel,  a  set  of  problems 
of  the  type  to  be  discussed  here  have  taken  on  greater  importance.   Since 
many  of  the  most  interesting  problems  in  this  field  have  system  equations 
for  which  the  methods  of  the  Calculus  of  Variations  fail,  new  mathemati- 
cal methods  are  required  for  their  solution. 

Suppose  that  the  system  to  be  investigated  may  be  described  by  an 
equation  of  the  form 


~=G(X,U)  ,      X(0)  =C  (1-D 


where  X(t)  is  the  state  variable  and  U(t)  is  the  control  variable.*  The 
problem  then  is  to  determine  U(t)  such  that  the  system  is  optimized 
according  to  some  specified  criterion  or  cost  function. 
Two  types  of  process  are  of  particular  importance: 

(1)  A  desire  for  X(t)  to  remain  as  close  to  a  prescribed  path  as 
possible  throughout  the  duration  of  the  process  and, 

(2)  A  desire  for  the  final  value  of  X(t)  to  be  a  prescribed  value 
or  in  a  prescribed  state  (Terminal  Control). 


rEqn.  (1-1)  is  written  in  one  dimension  and  this  investigation  will 
be  confined  to  systems  of  the  first  order.   However,  all  methods  and 
techniques  employed  may  easily  be  extended  to  n-th  order  systems  where 
X  and  U  would  be  replaced  by  the  vectors  X  and  U. 


The  type  of  problem  as  set  forth  in  case  (1)  is  the  type  investigated 
here  and  is  the  "Free  Right  End"  problem  as  described  by  L.  I.  Rozonoer 
in  Section  I  of  his  paper  dealing  with  Pontryagin's  Maximum  Principle  HI. 

A  non-linear  system  was  hypothesised  and  methods  of  obtaining  a 
control  policy  which  would  provide  an  optimum  system  according  to  some 
specified  criterion  were  investigated. 

The  system  chosen  for  the  investigation  may  be  described  by  a  non- 
linear differential  equation  of  the  first  order  which  has  the  following 
functional  form: 

X  =  f(X,U)   ,       X(0)  =  C  C1"2) 

in  which  one  may  consider  X  as  representing  the  error  or  deviation  from 
a  prescribed  path  and  U  as  representing  the  control  effort  applied.   The 
problem  was  to  determine  the  control  policy  U,  such  that  some  criterion 
function  of  the  form 

J  =  \  G(X,U)dt  (1-3) 

would  be  minimized.   Such  functions  as  described  by  Eqn.  (1-3)  with  the 
constraint  of  Eqn.  (1-2),  which  evaluate  the  performance  of  a  system  are 
commonly  referred  to  in  the  literature  as  "criterion"  or  "cost"  functions. 
Sometimes  they  can  be  expressed  simply  in  terms  of  integrals.  At  best, 
the  choice  of  a  cost  function  J,  is  a  compromise  between  a  desired  crite- 
rion of  goodness  for  the  control  design  and  one  which  leads  to  a  more 
tractable  mathematical  analysis. 

The  specific  system  investigated  is  described  by  the  first  order, 
non-linear  differential  equation 

X  =  AX  +  B  Sin  +  CU  (1-4) 


with  the  constants  having  the  following  values ; 

A  -  -0.1  B  =  +1.0  C  -  +.25 

Several  different  cost  functions,  subject  to  the  constraint  of 
Eqn.  (1-4)  were  proposed  for  the  investigation; 


=  (V+u2) 


dt  (1~5) 


(  Of 


-  (V>< 

'0 


J  =  \  (Tf)dt  (1-7) 

'o 

The  problem  of  obtaining  a  control  policy  to  effect  the  minimization 
of  these  criteria  (each  of  which  leads  to  a  different  set  of  adjoint 
equations)  will  be  treated  in  detail  below. 

Qualitatively,  minimizing  Eqn.  (1-5)  would  correspond  to  minimizing 
the  error  and  the  control  effort  to  the  system.  Minimizing  Eqn.  (1-6) 
would  correspond  to  minimizing  only  the  error,  while  minimizing  Eqn.  (1-7) 
would  minimize  only  the  control  effort  of  the  system  with  no  regard  for 
the  error. 

As  mentioned  above,  this  investigation  deals  with  the  "Free  Right 
End"  problem.  As  a  result,  the  time  interval  over  which  the  integrals, 
Eqs.  (1-5),  (1-6)  and  (1-7),  were  evaluated  was  fixed  at  10  seconds. 

By  the  very  nature  of  the  formulation  of  the  problem  (minimization 
of  the  integral   J  =  f  G(X,U)dt   with  the  constraint  X  =  f  (X,U),  having 

Jo 

X(0)  specified,  but  not  specifying  the  value  of  X  at  the  final  time  T)s  it 
falls  within  the  province  of  the  Calculus  of  Variations.   Occasionally 
such  classical  methods  can  be  employed  to  determine  the  optimal  control 
policies, 

The  authors  however,  attempted  to  apply  the  work  of  a  Russian  control 


theorist,  L.  S.  Pontryagin,  as  it  is  presented  in  a  series  of  articles 
by  L.  I.  Rozonoer  /1,2,3/.   With  the  aid  of  these  papers  and  the  CDC 
1604  digital  computer,  satisfactory  control  policies  were  obtained. 

Briefly  recapitulating,  the  purpose  of  this  thesis  is  to  use  the 
digital  computer,  employing  various  programming  methods,  to  search  for, 
compute,  and  design  an  optimal  control  which  will  satisfy  the  various 
criteria  (cost  functions)  as  applied  to  Eqn.  (1-4). 


2.   The  uncontrolled  system. 

The  uncontrolled  system  as  described  by  the  following  equation 


-0.1X  +  Sin  X 


(2-1) 


was  investigated  first.  A  knowledge  of  the  behavior  of  the  system  with- 
out restraint  or  control  applied  was  deemed  necessary  in  order  that 
investigations  of  the  system  with  control  could  be  properly  interpreted 
and  understood. 

The  Donner  Analog  computer  (Model  3100)  was  employed  in  order  to 
obtain  the  phase  plane  plot  of  velocity  (X)  versus  trajectory  (X).   Fig, 
(2-1)  shows  the  computer  diagram  used  to  obtain  the  phase  plane  desired. 


20  Sin(t) 


_J 


-Y  =  -2t  +  20  Sin(t) 


FIG.  (2-1).   Computer  diagram  utilized  in  order  to  obtain  the 
Phase  Plane  (X  vs.X)  for  the  uncontrolled  system. 

As  no  sine  function  generator  was  available,  the  sine  function  was 

simulated  as  shown  within  the  dotted  rectangle.* 


*The  sine  function  was  simulated  by  constructing  the  analog  of  the 
second  order  differential  equation  X  4  X  =  0  whose  solution  is  known  to 
be  X  =  sin(t). 


The  output 


-Y  =  20  Sin  t  -2t 


(2-2) 


was  recorded  on  an  X-Y  Plotter  utilizing  the  real  time  variable  input  to 
the  X  coordinate.   Fig.  (2-2)  shows  the  results  of  this  investigations, 
the  phase  plane  plot  for  the  uncontrolled  system. 


FIG.  (2-2).   Phase  Plane  Plot  of  the  uncontrolled  system  showing 
Velocity  (X)  on  the  ordinate  plotted  against  Trajectory  (X)  on 
the  abscissa. 

Since  the  computer  output  was  in  terms  of  Y  and  Times  it  was  neces 

sary  to  convert  these  coordinates  to  the  desired  coordinates  X  and  X. 

In  order  to  scale  the  abscissa  it  was  necessary  to  find  the  point  on 

the  curve  corresponding  to  rt  radians.   Since  the  curve  as  described  by 

Eqn.  (2-2)  is  a  linear  combination  of  a  sine  curve  and  a  straight  line 


U  =  -Kt 
6 


(2-3) 


the  points  where  the  resultant  curve  intersects  the  straight  line  are 
multiples  of  k.   The  points  on  the  abscissa  lying  directly  above  these 
points  of  intersection  are,  therefore,  also  multiples  of  rt  radians.   In 
this  manner  the  X  coordinate  was  transformed  from  a  time  scale  to  a 
distance  scale.   (Since  only  qualitative  results  as  to  velocity  were 
desired  the  Y  fxis  was  not  rescaled  and  no  transformation  of  coordinates 
was  made). 

Investigation  of  Fig,  (2-2)  yielded  the  following  interesting  informa- 
tion concerning  the  existance  of  equilibrium  points:  There  are  two  stable 
equilibrium  points  located  at  X  =  2.84+  and  at  X  -   8.41+  and  there  are  two 
unstable  equilibrium  points  located  at  X  =  0.0  and  at  X  =  7.02+.  Thus, 
the  uncontrolled  system  may  be  expected  to  move  toward  one  of  the  stable 
equilibrium  points  depending  on  the  initial  value  of  X  as  indicated  by 
the  arrows  on  the  curve  of  Fig.  (2-2). 

In  addition  to  the  above  investigation,  the  CDC  1604  digital  computer 
was  employed  as  an  additional  method  of  determining  the  equilibrium  points. 
By  choosing  selected  values  of  initial  X  (both  positive  and  negative)  and 
solving  the  differential  equation  (by  the  Runge-Kutta  method),  equilibrium 
points  were  obtained  which  agreed  favorably  with  those  obtained  by  the 
graphical  analysis  of  the  analog  solution.   (Fig.  (1-2)  and  Fig.  (1-3) 
in  Appendix  I  are  the  computer  graph  solutions  obtained  from  this  invest i= 
gation.   Tabular  output  also  was  obtained  but  is  not  included  here  due 
to  the  volume  obtained  and  due  to  the  relative  unimportance  of  this  sort 
of  detail). 

In  order  to  obtain  values  of  the  equilibrium  points  to  an  accuracy 
greater  than  that  possible  by  either  of  the  preceeding  two  methods, 
standard  mathematical  methods  were  utilized.  The  derivative  was  set 


equal  to  zero  and  the  values  of  the  equilibrium  points  were  obtained  to 
three  decimal  point  accuracy  with  the  aid  of  sine  tables,  a  desk  cal- 
culator and  iterative  proceedures.   The  results  obtained  again  corres- 
ponded to  and  verified  the  earlier  results,  yielding  stable  equilibrium 
points  at  X  =  2.852+  and  at  X  =  8.416+  with  two  unstable  equilibrium 
points  at  X  =  0.00  and  X  ■  7.068+.   Only  positive  equilibrium  points 
were  calculated  as  this  investigation  was  to  be  limited  to  positive 
values  of  X. 

Interpreting  the  significance  of  the  sum  of  these  investigations, 
it  is  possible  to  define  the  regions  of  stability  and  instability  for  the 
uncontrolled  system  as  shown  below  in  Fig.  (2=3) 


00 
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FIG.  (2-3)   Regions  of  stability  and  instability  for  the 
uncontrolled  system. 


3.   Maximum  principle  In  optimal  system  theory. 

As  stated  above,  in  an  effort  to  obtain  a  solution  to  the  problem 
of  obtaining  optimal  controls,  the  authors  utilized  the  "Maximum 
Principle"  as  hypothesized  in  1956  by  L.  S.  Pontryagin  on  the  basis  of 
work  performed  by  him,  V.  G.  Boltyanskii  and  R.  V.  Gamkrelidze.  /4/ 
Pontryagin' s  Maximum  Principle  is  presented  in  some  detail  by  L.  I. 
Rozonoer  in  a  series  of  articles  appearing  in  "Automatika  i  Telemakhanika' 
in  1959.   ("Automation  and  Remote  Control"  presents  an  English  transla- 
tion of  this  Russian  Journal.  /1.2.3/).   A  brief  resume  of  the  most 
important  concepts  will  be  made  here  in  order  that  a  common  ground  for 
further  discussion  may  be  establ  shed. 

Rozonoer' s  papers  deal  with  the  problem  of  obtaining  a  control  which 
is  optimum  in  some  sense  for  a  system  which  may  be  described  by  a  set  of 
differential  equations  of  the  n-th  order:* 

X.  =  f(X.,U.,t)     i=l,8.o,n  (3-1) 

where  X.,  ...,  X  are  the  parameters  of  the  system  and  U„,...,U  are 
the  positions  of  the  controlling  elements. 

Rozonoer  shows  in  his  paper  that  the  problem  of  optimizing  the 
system  with  respect  to  an  integral 


fc(X,U)dt 


(3-2) 
leads  to  the  problem  of  optimizing  with  respect  to  coordinates.**  At  this 


*Any  n-th  order  differential  equation  may  be  expressed  in  terms  of 
n  first  order  differential  equations  involving  n  variables. 

**L.  I.  Rozonoer,  L.  S.  Pontryagin' s  Maximum  Principle  in  the  Theory 
of  Optimum  Systems,  "Automation  and  Remote  Control",  Vol  20,  p!291. 


point  an  additional  variable  may  be   introduced 


fp(X1,...,Xn;U1,.98,Ur;t)dt     : 


Xn+l(t)=     K VU1 Vt)dt     :     Xn+1(0)=0  (>3) 


which  allows  one  more  differential  relation 


Wl'*1! Vui V*) 


(3-4) 


to  be  added  to  the  system  specified  in  Eqn„  (3-1).   The  problem  thus 
becomes  one  of  optimizing  the  n+lst  system  coordinate  at  the  final 
moment  of  time. 

Specifically,  the  problem  of  optimizing  a  linear  function  of  the 
final  values  of  all  the  coordinates  of  the  system,  that  is,  the  quantity 


s  =  2__  c 

k=l 


A(T)  (3-5) 


where  C,  are  certain  constants,  is  developed  in  detail  in  Rozonoer's 
paper.   In  the  discussion  of  this  problem,  the  theory  is  first  developed 
for  the  case  in  which  the  right  end  of  the  trajectory,  X(t),  is  not  fixed, 
That  is,  there  are  no  restrictions  imposed  on  the  final  values  of  the 
coordinates.* 

A  variable  vector,  P(t)  =  P.(t),...,  P  (t),  which  has  a  direction 
at  time  t  =  T  opposite  to  the  direction  of  the  vector  C  =  C.  ,.„.,,  C  is 
introduced  at  this  point 

P.(T)=-C.     i  =  l n  (3-6) 

It  is  assumed  that  the  moduli!  of  the  vectors  P(T)  and  C  are  equal. 


*For  the  entire  problem  considered  here,  the  final  coordinate  was 
left  free.   Only  the  duration  of  the  problem  was  fixed  at  10  seconds. 


10 


The  variables  P.(t')  are  subject  to  a  set  of  differential  equations 
n   5f  (X  ,...,X  ;IL,...,U  ;t) 

p.(t)=-7>  JLJ SLi L-^u-i b    (3-7) 

1      is        bx± 

(Let  us  note  here  that  if  any  control,  U(t),  is  given,  the  vector  P(t) 
is  uniquely  determined  from  Eqs.  (3-7),  where  conditions  Eqs.  (3-6) 
play  the  role  of  boundary  conditions.)   From  Eqs.  (3-5)  and  (3-6)  bound- 
ary conditions  for  the  final  values  of  the  adjoint  variables,  P„,  may  be 
obtained. 

Rozonoer  continues,  and  points  out  that  if  the  following  expression 
is  formed, 

H(X,P,U,t)  =  Yi   psxs  O-s) 

s=l 

that  Eqs.  (3-2)  and  (3-7)  may  be  written  in  the  form 

V%    •  *i--^..   '-1 ' 

Since  the  X.(0)  are  specified  in  the  problem  statement  and  the 
P.(T)  may  be  found  from  Eqs.  (3-6)  and  (3-5),  the  boundary  conditions 
on  Eqn.  (3-9)  are 


X.(0)=X°     ,  P.(T)=-C.     i  =  l,...,n       O-10) 

The  H-function  is  analogous  to  the  Hamiltonian  in  analytical 
mechanics,  and  the  vector  P(t)  to  the  impulse  vector.   Rozonoer  proves 
in  his  paper  that,  according  to  the  maximum  principle,  the  H-function 
must  be  maximum  (minimum)  in  U  for  any  values  of  X  and  P  in  order  to 
obtain  the  optimum  condition. 

Now  well  known  principles  of  Calculus  are  resorted  to  and  the 
partial  derivative  of  the  H-function  is  taken  with  respect  to  U  and  set 


11 


equal  to  zero, 

(3-11) 


fc-° 


Eqn.  (3-11)  yields  a  relation  between  the  control  U  and  the  adjoint 
variable  P. 

Since  the  cost  function  to  be  minimized  has  been  defined  as  a 
linear  combination  of  the  final  values  of  the  coordinates  of  the  object 
system,  Eqn.  (3-5),  and  this  function  added  to  the  system  as  the  n+lst 
coordinate,  it  may  be  seen  from  Eqn.  (3-9)  that 


P  .  =  0  (3-12) 

n-fl 

in  all  cases.  This  allows  immediate  integration  of  P  .  yielding 

n+i 

P  ,  =  CONSTANT  (3-13) 

n+l 

Since  P  ,  is  constant  and  we  know  from  Eqn.  (3-10)  that  P  , (T)  ■  -C  , 
n+l  H        v    '       n+lv  '  n+l 


we  see  immediately  that 


Also  note  that  C  -  will  always  be  equal  to  one  and  that 


Thus  we  see  that  P  1  =  -1  which  allows  some  simplification  in  Eqs.  (3-9) 
As  a  result  of  the  above  manipulations,  2n-l  differential  equations 
are  obtained,  the  solution  of  which  yield  the  desired  optimum  control. 
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2  2 
4.   The  cost  function  J  =   \  (X  +U  )  dt. 


4.0 


The  first  system  investigated  was  one  in  which  the  integral 

Idt  (4-1) 


=  (V+u2)< 


was  to  be  minimized.   The  system  to  be  controlled  was  described  by  the 
differential  equation 

X  =  AX  +  B  SinX  +  CU  (4-2) 

for  which  the  constants  A,  B,  and  C  had  the  same  values  as  specified  in 
Section  2;   namely,  -0.1,  1.0,  and  0.25  respectively.   The  time  interval 
T,  was  fixed  at  ten  seconds.   (Some  investigation  as  to  the  effect  of 
extending  the  time  interval  to  20  seconds  was  investigated  and  is  dis- 
cussed below.) 

Since  the  investigation  of  the  uncontrolled  system  showed  several 
points  of  stable  and  unstable  equilibrium  for  initial  values  of  X 
between  zero  and  10,  an  initial  value  of  X  equal  to  20  was  chosen  in 
order  to  include  the  effects  of  the  equilibrium  points. 

The  problem  was  approached  according  to  the  principles  as  set  forth 
in  Rozonoer's  paper,  /l/,  dealing  with  Pontryagin's  Maximum  Principle. 

In  order  to  apply  Pontryagin's  Maximum  Principle,  it  is  necessary 
that  the  system  be  described  by  an  equation  of  the  form  specified  by 
Eqn.  (3-1);  i.e.,  a  first  order  differential  equation,  either  vector  or 
scalar.   Since  the  specific  system  under  investigation,  Eqn.  (4-2),  is  of 
first  order,  no  additional  manipulation  was  required. 

Then,  in  accordance  with  Rozonoer's  statement  that  optimizing  with 
respect  to  an  integral  leads  to  the  problem  of  optimizing  with  respect  to 
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the  final  values  of  the  coordinates,  the  cost  function  was  added  as 
the  n+lst  coordinate 


'0 


^  =    \  u  +ir;dt  (4-3) 


The  system  equations  now  have  the  form 


X.  =  AX.  +  B  SinX1  +  CU  ,.    . . 

11        1  (4-4) 

i2  =  x2  +  u2 


Next,  forming  the  final  value  functional  S,  Eqn.  (3-5) 


=  7^.7.01!)=  C2X2(T)  =  f  (X^+U2)dt  (4-5) 


the  values  of  the  constants  C.  are  obtained, 

1  o1  =  o 


c2  =  i 


(4-6) 


and,  since  we  know  from  Eqn.  (3-6)  that  the  adjoint  variables  are  equal 

to  but  of  opposite  sense  to  the  C.  at  the  final  time  T,  we  now  have  the 

boundary  conditions 

PX(T)  =  0 

(4-7) 
P2(T)  =  -1 

Forming  the  H-Function,  Eqn.  (3-9), 

H(p,x,u,t)  =  yp.f.  =  Y?±\  (4-g) 

and  making  appropriate  substitutions  for  the  X.  the  H-function  becomes 

H  =  P-^AX-j^  +  B  SinXj^  +  CU)  +  ^Xl+   U2)  .  (4-9) 

From  Eqs.  (3-9)  and  (4-9)  it  is  now  possible  to  obtain  the  adjoint 
equations  below: 
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P2  =  -(A  +  B  CosX1)P1  -  2X^2 

(4-10) 

p2  =  o 

Applying  the  principles  of  Calculus  to  the  H-function  in  order  to 

obtain  the  minimum  with  respect  to  the  control  U,  a  relation  between 

the  adjoint  variables  Pi  and  the  control  variables  is  found. 

CP 
$=CP1  +  2UP2=0        U.-l  (4-11) 

Noting  that  the  equation  for  P.,  Eqn.  (4-10),  is  readily  integrable, 
it  is  seen  that  P2(t)  is  equal  to  a  constant.   Since  P2(T)  =  -1,  Eqn. 
(4-7),  it  follows  that  P?  must  equal  -1  at  all  times.   Knowing  the  value 
of  P2(t),  it  is  possible  to  simplify  Eqn.  (4-11). 

U  =  gi  (4-lla) 

Finally,  making  the  appropriate  substitutions  for  the  constants 
A,  B,  and  C  and  for  the  adjoint  variable  P~(t)  and  replacing  U  by  P,/8, 
the  following  set  of  differential  equations  are  obtained: 


X1  =  -.2X±   +  Sin  Xx  +  ,03125  P1  (4-12) 


X2  =  Xl  +  64 


(4-13) 


P  =  -(.1  +  Cos  X1)P1  +  2XX  (4-14) 


The  solution  of  the  above  set  of  equations  is  the  problem  remaining  at 
this  point. 

Numerical  methods  and  the  CDC  1604  digital  computer  were  utilized 
in  order  to  obtain  a  solution  to  this  problem.   Specifically,  the  Runge- 
Kutta  method  of  evaluating  first  order  differential  equations,  as  program- 
med for  the  CDC  1604,  was  employed.   (See  Appendix  V). 
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In  order  to  obtain  solutions  to  the  set  of  differential  equations 
by  numerical  methods,  initial  conditions  are  required.   Since  initial 
conditions  are  known  only  for  the  X  variables,  (X-(0)=20,  X2(0)=Q),  and 
not  for  the  adjoint  variable  P~,  the  immediate  problem  becomes  one  of 
searching  the  adjoint  space  for  the  initial  conditions  which  will  yield 
the  desired  solution.* 

The  first  method  proposed  in  an  effort  to  obtain  P*(0)  was  one  in 
which  several  P(0)  would  be  selected  in  an  effort  to  bracket  P*(0).   Then, 
through  some  selective  iterative  scheme,  the  bracket  size  would  be  reduced 
and,  eventually,  the  P*(0)  yielding  THE  OPTIMUM  system  could  be  obtained. 
(This  method  of  solution  is  commonly  called  a  "hill  climbing"  technique). 
There  was,  however,  no  evidence  that  the  variation  of  the  cost  function, 
J,  versus  the  initial  values  of  P  would  be  smooth,  and  the  possibility  of 
"homing  in"  on  some  local  minimum  rather  than  the  true  minimum  was  present, 
(The  investigation  of  the  system  showed  that  local  minima  did  exist  in 
many  cases.) 

In  view  of  the  above,  the  authors  decided  on  another  approach  to  the 
problem  of  obtaining  P*(0).  A  reasonable  range  of  values  for  P*(0)  was 
guessed.   Then,  the  differential  equations  for  the  system  were  evaluated 
for  Xx(0)  =  20,  X2(0)  =  0,  and  P^O)  ranging  from  -250  to  +50.  This 
investigation  yielded  values  of  cost  function  for  300  integral  values  of 
P(0)  as  shown  in  Fig.  (4-1). 

A  brief  comment  on  the  details  of  programming  this  problem  might  be 


*It  should  be  pointed  out  here  that  any  solution  generated  from  Eqs. 
(4-12),  (4-13)  and  (4-14)  will  be  optimal  with  respect  to  the  initial 
conditions  chosen.   However,  of  all  these  optimal  solutions,  there  is  one 
"best"  solution.   The  P(0)  which  yields  THE  OPTIMUM  solution  will  be 
designated  P*(0).   (The  asterisk  when  applied  to  any  other  variable  will 
likewise  refer  to  THE  OPTIMUM  system.) 
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FIG.  (4-1).   Cost  function  vs.  Initial  value  of  P. 
X(0)  =20  J  =  C(X2+U2)dt 

appropriate  here.   Since  it  had  been  decided  to  obtain  a  P*(0)  by 
calculating  the  cost  function  for  several  P(0),  the  next  question  to 
resolve  was  in  regard  to  the  range  of  values  of  P(0)  to  investigate.   It 
was  suggested  that  a  first  guess  might  be  that  the  initial  control,  0(0), 
should  be  at  least  as  large  as  X(0)  but  of  opposite  sign.   Utilizing  Eqn. 
(4-lla)  and  recalling  that  X(0)  ■  20,  a  first  guess  pucs  P(0)  at  about  -160. 
In  an  effort  to  allow  for  some  uncertainty  and  also  to  obtain  information 
as  to  the  behavior  of  the  cost  function  for  P(0)  not  in  the  vicinity  of 
P*(0),  a  range  of  F(0)  was  selected  from  -250  to  +50  as  stated  above. 
Fig.  (4-1)  was  the  result  of  this  initial  investigation.   (Not  shown  in 
Fig.  (4-1),  but  obtained  in  a  tabular  output  were  the  final  values  of  the 
trajectory,  X(T).) 

Locating  the  P*(0)  which  yields  the  minimum  cost  function  and  then 
repeating  the  above  investigation  in  the  neighborhood  of  the  first  deter- 
mined P*(0)  +€  allows  P*(0)  to  be  determined  to  any  degree  of  accuracy 

desired. 
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The  major  disadvantage  of  the  above  method  is  the  time  required  to 
obtain  P*(0)  to  a  reasonable  accuracy.  In  excess  of  11  minutes  was  re- 
quired for  the  above  investigation.  (It  was  not  possible  to  obtain  tra- 
jectories or  control  policies  during  this  first  investigation  due  to 
storage  limitations  of  the  computer.)  At  that,  the  resulting  P*(0)  was 
accurate  only  to  an  integral  value. 
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FIG.  (4-2).   Optimal  control  policies,  U(t),  showing  variation 
with  the  accuracy  to  which  P*(0)  is  computed  for  J  ■  J( 


(X2+U2)dt, 


It  should  be  noted  that  the  accuracy  to  which  P*(0)  was  determined 
affected  both  the  control  policy  and  the  trajectory.   Fig.  (4-2)  shows 
the  effect  on  control  that  the  accuracy  of  the  P*(0)  had. 

The  OPTIMAL  Control  Policy  obtained  from  the  above  computations 
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yielded  a  cost  function  of  1099.52  versus  a  cost  function  of  1978,09  for 
the  uncontrolled  system,  an  improvement* of  44.8%. 
20 


FIG.   (4-3).  Trajectory  vs.  Time  for  the  OPTIMAL  system  (solid 
curve)  and  the  uncontrolled  system  (dotted  curve). 

Fig.  (4-3)  shows  the  trajectories  obtained  for  the  uncontrolled 

system  and  the  optimally  controlled  system.   Note  that  the  trajectory 

obtained  for  the  optimal  case  as  shown  in  Fig.  (4-3)  appears  to  approach 

a  constant  value  of  about  2.7  rather  than  a  final  value  of  zero.   However; 

recall  that  the  cost  function  to  be  minimized  here  is  a  combination  of 

error,  X,  and  control  effort,  U,  and  also  that  the  uncontrolled  system 

had  a  stable  equilibrium  point  at  X  =  2.84.   One  might  assume,  therefore, 


*Percent  improvement  is  defined  as 

Percent  improvement  =  J  (uncontrolled)  -  J  (controlled) 

J  (uncontrolled) 
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that  the  effort  necessary  to  reduce  the  error  below  this  value  would 


exceed  any  reduction  in  the  cost  function  that  would  result  from  such  a 

TIME  (Sec) 
4  6  8  10 


reduction  in  X. 

2 
4 


FIG.   (4-4).   Optimal  control  policy  for  X(0)  =  20  showing 
local  maxima  and  minima. 

Fig.  (4-4)  shows  the  optimal  control  policy  obtained  for  this  system 
with  X(0)  =  20.   It  is  very  interesting  to  note  that  the  points  of  local 
maxima  and  minima  occur  when  the  trajectory,  Fig.  (4-3),  has  a  value  such 
that  the  velocity,  X,  of  the  system  is  also  very  close  to  a  maximum  or 
minimum.   (See  Fig.  (2-2)  which  is  the  velocity-displacement  phase  plane 
for  the  uncontrolled  system). 

Having  investigated  the  system  with  the  initial  value  of  X  set  at 
20,  the  next  part  of  the  investigation  involved  the  determination  of  the 
effect  of  the  initial  value  of  X  on  the  P*(0),  the  optimum  control  policies 
and  the  optimum  trajectories.   Consequently,  the  systems  having  initial 
values  of  X  equal  to  15,  10,  and  5  were  investigated.   The  P*(0)  in  each 
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case  was  obtained  in  a  similar  manner  to  that  described  for  the  system 
in  which  X(0)  =  20.   Optimal  controls  and  optimal  trajectories  were 
obtained  in  each  case. 

Upon  investigation  of  the  resulting  optimal  control  policies,  it 
was  noted  that  they  were  remarkably  similar  in  appearance'  to  the  optimal 
control  policy  for  the  system  having  an  initial  X  equal  to  20.   Fig.  (4-5) 
shows  the  control  policies  for  the  four  systems  (X(0)  =  20,  15,  10,  5). 
Note  that  the  correspondence  is  very  good  for  X(0)  -  10  and  15  until  the 
latter  part  of  the  common  time  interval  at  which  point  the  three  policies 
begin  to  separate.   Of  course, this  correspondence  only  holds  if  the  origin 
of  the  time  axis  for  the  systems  having  initial  X  equal  to  10  and  15  are 
moved  to  the  right  as  indicated.   Notice,  however,  that  the  control  policy 
for  the  system  having  an  initial  X  equal  to  5  does  not  seem  to  correspond 
at  all. 

In  spite  of  this  latter  result,  the  question  arose  as  to  whether  the 
optimal  control  policy  for  the  X(0)  =  20  system,  in  conjunction  with  its 
trajectory,  might  provide  the  P*(0)  for  all  initial  values  of  X  between 
20  and  the  lower  limit  of  the  trajectory  of  about  2.7. 

Consequently,  systems  having  initial  values  of  X  of  17.5,  12.5,  7.5, 
2,5,  and  1.0  were  investigated  for  P*(Q).   Having  obtained  the  P*(0)  for 
these  systems,  the  calculated  P*(0)  were  ploted  on  a  curve  of  optimal 
control  versus  optimal  trajectory  for  the  system  having  initial  X  equal 
to  20,  Fig.  (4-6).   One  immediately  notices  the  close  correspondence 
except  for  those  system  having  initial  X  equal  to  5,  2.5,  and  1.0.   Since 
the  optimal  trajectory  for  the  X(0)  =  20  system  never  got  below  a  value 
of  about  2.7,  it  is  reasonable  to  assume  that  it  would  not  be  possible  to 
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FIG.  (4-5)  Control  vs.  Time  for  various  initial  X  showing  corres- 
pondence except  during  the  latter  port  of  the  common  time  interval  and  showing 
the  non-correspondence  for  initial  X  equal  to  5.0. 
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PIG. (4-6).  Optimal  Control  plotted  against  Optimal  Trajectory  for  the 
system  having  initial  X  equal  to  20.  Indicated  points  are  the  P*(o)  obtained 
for  systems  having  inital  X  equal  to  the  values  indicated.  For  T  =  10  sees. 
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obtain  a  P*(0)  for  initial  X  which  were  less  than  this  value,  i.e.,  not 
on  the  curve. 

As  a  result  of  the  above  investigation,  the  possibility  was  suggested 
that  if  the  system  having  initial  X  equal  to  20  were  allowed  to  run  for  20 
seconds  and  a  new  optimal  control  policy  and  trajectory  obtained  for  this 
system,  that  it  might  be  that  the  control  policy  for  the  X(0)  =  5  system 
would  match  and  that  this  new  U*  vs,  X*  would  provide  the  P*(0)  for  the 
values  of  X  less  than  2.7. 

The  X(0)  =  20  system  was  allowed  to  run  for  20  seconds  and  a  new 
optimum  control  and  trajectory  were  obtained  as  well  as  a  plot  of  U*  vs. 
X*.  When  the  P*(0)  for  the  system  investigated  above  were  plotted  on  the 
new  U*  vs  X*  curve,  Fig.  (4-7),  very  close  correspondence  was  noted  in  all 
cases  including  the  ones  for  initial  X  of  5,  2.5  and  1.0.   One  must  keep 
in  mind  that  by  extending  the  time  interval  to  20  seconds  the  problem 
became  a  new  problem  entirely.   Even  so,  the  results  of  this  investiga- 
tion show  that  the  U*  vs.  X*  curve  for  X(0)  =  20  over  the  extended  time 
interval  may  be  used  to  obtain  an  initial  guess  for  the  P*(0)  for  values 
of  X(0)  between  20  and  zero.   Modification  of  the  initial  guess  in  order 
to  obtain  the  true  P*(0)  would  be  a  much  easier  problem  than  the  problem 
of  obtaining  P*(0)  from  scratch. 

The  results,  as  stated  above  and  as  shown  in  the  various  figures, 
prove  conclusively  that  it  is  possible  to  obtain  numerical  solutions  to 
the  problem  of  obtaining  an  optimal  solution  to  the  non- linear  system  by 
means  of  Pontryagin's  Maximum  Principle. 

It  is  to  be  pointed  out  that  a  very  distressing  feature  of  this 
method  of  solution  is  the  excessive  amount  of  time  involved  in  obtaining 
a  satisfactory  solution  to  the  problem  compared  v ith  the  problem  time 
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FIG. (4-7) •  Optimal  Control  plotted  against  Optimal  Trajectory  for  the 
system  having  initial  X  equal  to  20.  The  time  interval  was  extended  in 
this  case  to  20  seconds.  Indicated  points  are  the  P*(o)  obtained  for  the 
systems  having  initial  X  equal  to  the  values  indicated. 
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of  10  or  20  seconds.   (This  would  siggest  some  investigation  into  the 
possibility  that  more  sophisticated  programming  of  the  CDC  1604  digital 
computer  might  result  in  a  shorter  solution  time  on  the  computer.) 

Finally,  one  additional  fact  should  be  borne  in  mind.   Although  the 
digital  computer  was  utilized  successfully  to  obtain  the  optimal  control 
function  (specified  at  intervals  of  0.1  seconds),  building  a  physical 
component  to  duplicate  this  control  is  quite  another  matter.   Some  com- 
promise would  probably  be  necessary.   Perhaps  a  square  pulse  of  U  (as  in 
the  output  of  a  zero  order  hold  in  a  sample  data  system)  corresponding  in 
width  or  modulus  to  the  peaks  shown  in  Fig.  (4-5)  would  suffice  for  engine- 
ering purposes.   Another  possibility  would  be  to  feed  some  percentage  of 
the  value  of  the  trajectory  back  to  the  controller  (negative  feedback). 

In  any  case,  each  type  of  control  would  have  to  be  evaluated  and 
the  cost  function  compared  to  the  optimum  one.   It  then  becomes  a 
question  of  deciding  whether  or  not  the  results  are  sufficiently  "opti- 
mum" for  the  engineering  purpose  in  mind.   Perhaps  the  only  utility  of  an 
investigation  such  as  this  is  in  providing  the  "ideal"  with  which  the 
engineer  may  judge  the  performance  of  the  system  he  has  designed. 

(Additional  graphs  of  the  optimal  control  policies  and  the  optimal 
trajectories  for  the  above  mentioned  initial  X  supplementing  and  support- 
ing the  above  data  and  conclusions  may  be  found  in  Appendicies  I  and  II.) 


26 


5.      The   cost    function    J  =\    U  dt. 


£ 


■ !?  «• 


The  second  system  investigated  was  one  in  which  the  integral 

(5-1) 
was  to  be  minimized  with  the  same  object  equation  as  in  Section  4. 

X  =  AX  +  B  Sin  X  +  CU  (5-2) 

Proceeding  in  the  same  manner  as  in  Section  4,  introduce  the  vari- 

ables  x1  =  X(t) 


£' 


*2  "  h  *"   dt  ' 
Substituting  these  variables  in  to  Eqs.  (5-1)  and  (5-2),  and  using 

the  differential  form  yields  the  system  equations 

Xx  =  AX1  +  B  Sin  X1  +  CU  (5-3) 

X2  =  U2  .  (5-4) 

By  forming  the  final  value  functional 

S  =  X2(T)  =  y"0^±   =  C^  +  C2X2  (5-5) 

the  C  vector  is  obtained.   In  order  to  minimize  the  final  state  of  the 
adjoint  equations,  recall  that  P.(T)  must  equal  to  -  C..   Thus  the 
following  boundary  conditions  for  the  adjoint  space  are  obtained 

P1(T)  =  -C1  =  0 

P2(T)  -  -C2  -  -1  (5"6) 

Forming  the  H- function  as  before 

H  =  P1(AX1  +  B  Sin  Z±  +  CU)  +  ?£?  (5-7) 
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and  performing  the  partial  differentiations  as  defined  by  Eqs.  (3-9), 

the  following  set  of  simultaneous  equations  is  obtained, 
X  =  AX  +  B  Sin  Xx  +  CU 

X  -u2 

A2  "  U  (5-8) 


Pj^  =  -(A  +  B  Cos  X1)P1 


p2  =  o 

By  differentiating  H  with  respect  to  U  and  setting  equal  to  zero, 
a  relation  between  the  control  U,  and  the  adjoint  variables  is  obtained^ 

CP 

Integrating  the  last  of  Eqs.  (5-8),  one  is  able  to  evaluate  P«(t) 
as  before. 


And  since 


then 


P_  =  Constant. 


P2(T)  =  -  1 


P2(t)  =  -  1 


Now  substituting  values  for  P_  and  C  into  Eqn.  (5-9),  the  relation 
between  P  and  U  becomes 

u=-i  (5-9a) 

o 

Now  Eqs.  (5-8)  may  be  arranged  in  a  form  suitable  for  solution 


Xx  =  -.U^  +  Sin  Xx  +  .03125  ^  (5-10) 

i,  =  P2  /  64  (5-11) 


Px  =  (.lXj  -  Cos  X1)P1  (5-12) 
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having  the  following  boundary  conditions: 

X1(0)  =  20  PX(T)  =  0 


X2(0)  =  0  P2(T)  -  -1 


(5-13) 


As  in  Section  4,  two  of  the  above  equations^,  Eqs,  (5-10)  and  (5-12), 
must  be  solved  simultaneously  after  finding  the  P*{Q)  which  minimizes  the 
cost  function,  Eqn.  (5-1).   But,  before  stating  the  results  of  this  in- 


f/ 


vestigation,  recall  that  the  integral  to  be  minimized  is  1  IT  dt  o 

'0 
Even  without  Pontryagin's  Maximum  Principle,  the  solution  is  immediately 

obvious.   U  must  equal  zero  if  the  cost  function  is  to  be  minimized,  since 
the  control  U,  is  independent  of  X  in  the  object  equation,  Eqn.  (5-1). 
Thus,  one  might  hope  for  a  reliable  check  on  the  answers  obtained  by  the 
methods  utilized  in  this  paper  in  the  solution  to  this  system. 

Proceeding  as  in  Section  4,  values  of  P(0)  between  -  250  and  +  50 
were  investigated  and  the  corresponding  cost  functions  computed.   The 
results  verified  the  expected  minimum  of  zero  for  P*(0)  =  0.0.   Fig.  (5-1) 
shows  the  minimum  cost  function  for  an  initial  P  equal  to  zero  which 
resulted  in  a  zero  control  for  the  entire  time  interval;  i.e.  the  un- 
controlled system.   As  was  stated  above,  this  is  a  reasonable  and  expect- 
ed solution. 


-£* 


Utilizing  the  criteria  J  =  \  u  dt  provided  a  check  on  the  solutions 

b 

obtained  in  synthesizing  an  optimal  system  by  the  maximum  principle. 
The  results  obtained  from  the  solution  of  Eqs.  (5-10) >  (5-11)  and  (5-12) 
verified  the  predictions  made  through  inspection  of  the  H- function  and 
its  derivative,  Eqs.  (5-7)  and  (5-9). 

It  appears  from  the  above  discussion^  that  this  is  not  a  very  sensible 
way  to  use  UT  dt   as  a  cost  function.   A  more  meaningful  use  would  be  to 
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FIG.  .(5-1).   Cost  Function  vs.  Initial  Value  of  P.   X(0)  =  20. 
J  =   f(U2)dt,  showing  a  minimum  cost  function  J  s  0,  at  P(0)  =  0, 

apply  the  above  cost  function  in  a  problem  where  the  state  of  the  tra- 
jectory at  time  (T)  =  10  sees,  is  fixed,  say  X(T)  =  0.0. 

That  is,  to  even  have  a  reasonable  problem  we  must  have  a  fixed 
right  end  trajectory  and  also  fix  time,  (T) . 
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6.      The   cost    function   J  =  \  &     dt. 

Jo 

The  third  system  to  be  investigated  was  one  in  which  the  integral 

J  =  f  X2  dt  (6-1) 
was  to  be  minimized,  again  with  the  object  equation 

X  =  AX  +  BSinl  +  CU  (6-2) 

Proceeding  as  before,  it  is  seen  that  the  system  equations  are 

X_  =  AX,  +  B  Sin  X.  +  CU  ,,    _N 

11         1  (6-3) 

X  -X2 

*2  _X1 

Forming  the  H-function 

H(P,X,U,t)  =  P1(AX1  +  B  Sin  Xx  +  CU)  +  P^        (6-4) 

and  taking  suitable  differentials,  Eqs.  (3-9),  it  is  possible  to  obtain 

the  adjoint  equations.   Four  simultaneous  equations  are  thus  obtained 

which,  when  solved,  provide  the  solution  to  the  problem  of  obtaining 

the  optimal  system. 

X  =  AX  +  B  Sin  X  +  CU 

X  -X2 

A2  "  Al  (6-5) 


Pj^  =  -  P  (A  +  B  Cos  X  )  +  2P2X1 
P2  =  0 
From  the  final  value  functional,  Eqn.  (3-5), 

S  =£[  Cixi  "  C1X1  +  C2X2  =  )  A   dt         (6"6) 

the  C  vector  is  obtained.   As  before,  by  setting  P„(T)  equal  to  -C „ , 
the  terminal  boundary  conditions  for  the  P.(t)  are  obtained. 
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P1(T)  =  0 


P2(T)  =  -1 


(6-7) 


Since  the  last  of  Eqs.  (6-5)  may  be  integrated;,  it  is  seen  that 

P  (t)  =  CONSTANT  (6-8) 


and.  since 


it  is  seen  that 


P2(T)  =  -1 


P2(t)  =  -1  (6-9) 

Substituting  this  value  for  P?  into  the  H-  f  jinction,  and  then  dif- 
ferentiating with  respect  to  U  and  setting  the  resultant  relation  equal 
to  zero  in  order  to  find  the  minimum,  results  in 

P-jC  =  0  (6-10) 

Since  C  is  a  constant  equal  to  .25,  the  obvious  conclusion  must  be  that 

Px  =  0  (6-11) 

At  first  the  last  results  seem  quite  disheartening  as  not  relation 
between  U  and  the  adjoint  space  exist.   However,  this  might  have  been 
anticipated  since  the  H- function  is  linear  in  U  and,  as  such,  could  have 
a  minimum  only  at  one  of  the  boundaries. 

Reconsidering  the  Principle  of  the  Maximum  in  view  of  the  above 
results,  one  realizes  that  in  order  to  minimize  the  H  f anction,  Eqn. 
(6-4),  for  any  X  and  P,  U  must  be  as  large  as  possible  and  of  opposite 
sign  to  X  (since  the  constant  C  is  positive  and  it  is  desired  to  make 
that  term  negative).   But,  since  most  practical  systems  have  an  upper 
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maximum  on  the  amount  of  control  available,  U  would  then  have  to  assume 
the  largest  value  possible  under  such  a  constraint.*  In  mathematical 
notation,  the  above  result  would  be  expressed  as 

u  =  -senXlulmax  <6"12> 

As  a  result  of  the  above  developments,  the  only  system  equation 
needed  to  obtain  the  optimal  control  for  this  system  is 

X1  =  AX1  +  BSinX1-SgnX1  |ll|max  (6-13) 

Eqn.  (6-13)  was  solved  numerically  on  the  digital  computer  by  means 
of  the  Runge-Kutta  method  for  evaluation  of  Differential  Equations.   A 
computer  solution  was  desired  in  order  to  verify  the  conclusions  made 
on  the  basis  of  the  above  results  and  knowledge  of  the  uncontrolled  systems 
namely;  (1)  The  cost  function  should  decrease  with  increasing  values  of 
Control  U,  and  (2),   Once  the  system  approached  and  was  driven  through 
zero,  a  chatter  mode  should  be  obtained  for  the  remainder  of  the  interval. 

The  initial  computer  solutions  immediately  verified  that  the  Cost 
Function  did  decrease  as  the  absolute  magnitude  of  the  Control  was  in- 
creased.  However,  the  chatter  mode  did  not  appear  in  the  form  expected. 
Fig.  (6-1)  shows  that  for  a  Control  having  a  magnitude  of  25,  the  Tra- 
jectory approaches  zero  in  about  2.8  seconds  and  then  begins  to  increase 
for  a  short  while  and  then  is  driven  back  towards  the  zero  point.   This 
variation  of  the  trajectory  with  time  was  not  the  chatter  mode  expected. 


*This  result  is  similar  to  the  well  known  Bang-Bang  principle. 
Using  Pontryagin's  theory,  this  would  be  equivalent  to  minimizing  fdt 
used  as  the  final  value  functional. 
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FIG.  (6-1).   Trajectory,  X(t),  and  Control,  U(t)  plotted  as  a  func- 
tion of  time,  showing  the  chatter  mode  initially  obtained  from  the 
solution  of  the  Ja  dt   system. 

Note  also  that  the  control  remained  constant  at  -25  for  the  entire  inter- 
val.   The  authors  had  expected  that  the  negative  control  would  drive  the 
trajectory  past  the  zero  point,  change  sign  as  X(t)  became  negative  ac- 
cording to  Eqn.  (6-12),  and  then  drive  X(t)  positive..   The  cycle  should 
repeat  itself  until  the  end  of  the  interval.   These  results  were  not 
obtained  as  may  be  seen  in  Fig.  (6-1). 

It  was  suggested  that  perhaps  the  sampling  rate  (0.1  sec)  might  be 
too  coarse  to  show  the  chatter  and  that  "noise"  was  interfering  with  the 
computations.   Consequently,  a  sampling  rate  of  .025  seconds  was  tried9 
but  with  no  success.   Results  similar  to  those  shown  in  Fig.  (6-1)  were 
obtained. 

At  this  point,  the  program  being  utilized  in  the  solution  of  this 
system  came  under  close  scrutiny.   It  was  discovered  that  the  method 
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used  to  evaluate  the  differential  equation  (Runge-Kutta)  utilized  the 
last  X(t)  computed  in  the  interval  dt  to  set  the  sign  of  U(t)  for  the 
next  evaluation  of  X(t)  in  the  same  interval.   The  program  should  have 
been  using  the  sign  of  the  X(t)  computed  during  the  n-lst  time  interval 
to  control  the  sign  of  U(t)  during  the  nth  time  interval.   Because  of 
the  nature  of  the  Runge-Kutta  routine,  four  values  of  X(t)  were  being 
averaged  in  each  time  interval,  some  being  positive  values  and  some 
possibly  negative,  resulting  in  the  erroneous  results. 

The  program  was  revised  to  correct  the  above  deficiency  and  the 
results  obtained  showed  the  chatter  mode  and  the  alternating  nature  of 
the  control  as  expected.   See  Fig.  (6-2). 


TIME 


TIME 


-25 


FIG.  (6-2).   Sketch  showing  the  chatter  mode  for  X(t)  and  the 
alternating  characteristics  of  U(t)  for  the  system  having  jX  dt 
as  a  cost  function. 

The  results  obtained  using  the  modified  program  confirmed  the  earlier 
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results  that  the  cost  function  decreased  with  increasing  magnitude  of  U. 
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FIG.  (6-3).  A  curve  of  Cost  Function  plotted  against  |uj  showing 
the  decrease  in  cost  function  (J)  with  the  increase  in  |U|  . 

Through  investigation  of  this  system  employing  various  values  of  U, 
it  was  found  that  the  magnitude  of  the  control  effected  both  the  magni- 
tude and  the  frequency  of  the  chatter.   Additional  graphical  results 
may  be  found  in  Appendix  IV.   From  an  engineering  point  of  view,  it 
might  be  desirable  to  include  a  dead  zone  around  zero  in  order  to 
eliminate  or  reduce  the  chatter. 
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7.   Conclusions. 

In  this  investigation,  the  authors  have  inquired  into  and  shown  the 
feasibility  of  obtaining  a  numerical  model  for  the  optimal  control  of  a 
non-linear  system  described  by  a  first  degree  differential  equation  by 
means  of  the  theory  and  procedures  of  Pontryagin's  Maximum  Principle  as 
set  forth  by  L.  I.  Rozonoer  /l,  2,  T/.   The  investigation  was  limited  to 
the  "free  right  end"  problem  and  to  the  single  dimensional  case. 

The  methods  and  procedures  utilized  may  easily  be  extended  to  the 
n-dimensional  case  without  any  basic  change  in  theory  or  procedure. 

Various  cost  functions  have  been  investigated  and  some  conclusions 
verified  by  numerical  solution  utilizing  the  CDC  1604  digital  computer. 

Future  investigation  might  well  be  in  the  area  of  second  or  third 
order  systems  and  also  in  the  area  of  computer  programming  techniques. 
The  methods  and  programs  utilized  to  obtain  solutions  were  effective, 
but  not  very  efficient. 

The  theories  utilized  here  do  yield  solutions.   However,  other 
methods  (Dynamic  Programming  for  instance)  might  well  be  more  effective 
and  useful  to  the  investigator,  depending  upon  his  specific  needs. 
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APPENDIX  I 
GRAPHS  FOR  THE  UNCONTROLLED  SYSTEM 
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?IG.  (i-l)  Phase  plane  plot  showing  velocity  as  the  ordinate  against 
distance  as  the  abscissa  for  the  uncontrolled  system  having  initial  X  equal 
to  20.  This  figure  shows  the  graphical  transformation  of  the  abscissa  from 
real  time  to  distance.  (Scale i  i-,05r/Hxi&i     X-,2%   rad/line.) 
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FIG.  (1-2)  Uncontrolled  System  Trajectory  (X)  vs  Time,  for  11  positive 
values  of  X(o)  showing  the  two  stable  equilibrium  points,  X  =  2.85,  and  842 
for  the  positive  system. 
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FIG.  (l-3)  Uncontrolled  System  Trajectory  (x)  vs  Time,  for  11  negative 
values  of  X(0)  showing  the  two  stable  equilibrium  points,  X  =-2.85,  and  -8.52 
for  the  negative  system. 
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APPENDIX  II 
GRAPHS  FOR  THE  COST  FUNCTION  J  = 


(T  (x2+u2)dt 
Jo 
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GRAPHS  FOR  THE  COST  FUNCTION  J  = 


f(X2  4  U2)  dt 

Jo 
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FIG.  (ll-l)  Cost  Function,  J  =  )j,Xz+  U2)dt,  plotted  as  a  function  of 
initial  value  of  the  adjoint  variable  P,  showing  the  true  minimum  and  two 
local  minima  for  the  system  having  an  initial  X  =  20. 
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FIG.  (II-2)  Cost  Function,  J  =/(Xi+  U2')dt,  plotted  as  a  function  of 
initial  value  of  the  adjoint  variable  P,  showing  the  true  minimum  and  one 
local  minimum  for  the  system  having  an  initial  X  =  17.5. 
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FIG.   (II-3)  Cost  Function,  J  =/(Xa'+  Uz)dt,  plotted  as  a  function  of 
initial  value  of  the  adjoint  variable  P,  showing  the  true  minimum  and  one 
local  minimum  for  the  system  having  an  initial  X  =  15. 
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PIG.  (II-4)  Cost  Function,  J  =/(X2+  U2)dt,  plotted  as  a  function  of 
initial  value  of  the  adjoint  variable*?,  showing  the  true  minimum  and  two 
local  minima  for  the  system  having  an  initial  X  =  7.5. 
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PIG.  (II-6)  Cost  Function,  J  =Jx2+  Uz)dt,  plotted  as  a  function  of 
initial  value  of  the  adjoint  variable  P,  showing  the  true  minimum  and  two 
local  minima  for  the  system  having  an  initial  X  =  2,5. 
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PIG.  (II-7)  Cost  Function,  J  =j[(Xz+  U*)dt,  plotted  as  a  function  of 
initial  value  of  the  adjoint  variable  P,  showing  the  true  minimum  and  two 
local  minima  for  the  system  having  an  initial  X  =  1.0. 
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PIG.  (lI-8)  Optimal  Trajectory  (x)  plotted  against  Time,  for  two 
values  of  P*(0)  showing  the  variation  in  trajectory  as  a  function  of  the 
accuracy  to  which  P*(o)  is  determined;  for  the  system  having  the  Cost 
Function  J  =  PUS  U*)dt  and  initial  X  =  20. 
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PIG.  (II-9)  Optimal  Control  (U)  plotted  against  Time,  for  two  values 
of  P*(o)  showing  the  variation  in  control  as  a  function  of  the  accuracy-ito 
which  P*(o)  is  determined;  for  the  system  having  the  Cost  Function  J  =j(Xz+  U2-)^ 
And  initial  X  =  ?fL  rt 


and  initial  X  =  20. 
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FIG.  (II-IO)  Curves  of  Control  (u)  and  Trajectory  (X),  plotted 
against  Time,  for  the  optimal  system  having  the  Cost  Function  J  =  j7x2+  Ui)dt 
and  initial  X  =  15.  ° 
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PIG.  (il-ll)  Curves  of  Control  (u)  and  Trajectory  (x),  plotted 
against  Time,  for  the  optimal  system  having  the  Cost  Function  J  =  J  {xz+  U^dt 
and  initial  X  =  10. 
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FIG.  (II-12)  Curves  of  Control  (u)  and  Trajectory  (x).  plotted 
against  Time,  for  the  optimal  system  having  the  Cost  Function  J  =j[Xl+  U2)dt 
and  initial  X  =  5.  ' 
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FIG.  (il-  13)  Optimal  Trajectory  (X)  plotted  against  Time,  for  the 
system  were  T  final  was  extended  to  20  seconds,  having  the  Cost  Function 
J  =  /(X*  +  Ui)dt  and  initial  X  =  20.0.   (Note  when  the  time  was  extended 
from  10  to  20  seconds  that  the  trajectory  approached  zero  and  was  driven 
through  the  lower  stable  equilibrium  point,  X  =  2.85  by  the  optimum  controller,) 
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FIG.  (11-14)  Optimal  Control  (u)  plotted  against  Time,  for  the 
system  were  T  final  was  extended  to  20  seconds,  having  the  Cost  Function 
J  =4^2+  Tj2")dt  and  initial  X  =  20, 


57 


appendix  iii 

;t  2 
graphs  for  the  cost  function  j  =  i   u  dt 
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PIG.  (ill-l)  Cost  Function,  J  =  J0(U*)dt,  plotted  as  a  function  of 
initial  value  of  the  adjoint  variable  P,  showing  the  true  minimum  at 
?(0)  =  0.0  with  J  =  0,0  for  the  system  having  an  initial  X  =  20.0. 


59 


APPENDIX   IV 
GRAPH  X     dt 
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FIG.   (17-1)  Curves  of  Control  (u)  and  Trajectory  (x),  plotted  against 
Time,  for  the  optimal  system  having  the  Cost  Function  J  =j(Xz)   dt 
initial  X  =  20,  and  |U|MAX  =  5.  (Note  the  absence  of  a  chatter  mode,  as 
|UU**  was  not  sufficient  to  drive  the  trajectory  to  zero  in  the  alloted  time.) 
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FIG.  (IV-2)  Optimal  Trajectory  (x)  plotted  against  Time,  for  the 
system  having  the  Cost  Function  J  =  J  (Xz)dt  and  initial  X  =  20,  (Note 
establishment  of  a  chatter  mode.)  |uj^K  =  10, 


the 
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FIG.(lV-3)  Optimal  Control  (u)  plotted  against  Time, 
baring  the  Cost  Function  J  = j  (Z^)dt  and  initial  X  =  20.  (N 
ment  of  a  chatter  mode.)  !U|m^  =  10. 


for  the  system 
ote  the  establish- 
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PIG.  (IV-4)  Optimal  Trajectory^)  plotted  against  Time,  for  the 
system  having  the  Cost  Function  J  =  X(X2)dt  and  initial  X  =  20.   (Note  the 
establishment  of  a  chatter  mode.)  |U|*4x  =  25. 
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^  \s^^M§mmmm^^mmk^ 


FIG.  (TV-5)  Optimal  Control  (TT)   plotted  against  Time,  for  the  system 
having  the  Cost  Function  .T  SI  l*")dt  and  Lnitial  X  =  20.  (Note  .the  establish 
nent  cf  a  chatter  mode.)  {1Tlm<u  -   25* 
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FIG.  (lV-6)  Optimal  Trajectory  (x)  plotted  against  Time,  for  the 
system  having  the  Cost  Function  J  =  -UXz)dt  and  initial  X  =  20.  (Note  the 
establishment  of  a  chatter  mode.)  |Uj^A<=  50. 
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■7)  Optimal  Control  (u)  plotted  against  Time,  for  the  systei 
t  Function  J  =KX  )  ;  and  initial  X  =  20.  (Note  the  establj 
tor  mode.;      |UjWftx ..  50. 
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APPENDIX  Va 
SUBROUTINE  RUNGE - KUTTA 
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V.a.   The  Kunge-Rutta  method  of  evaluating  differential  equations. 

There  are  many  variants  of  the  Runge-Kutta  method,  but  the  most 
widely  used  one  is  the  following:   given 

y'  =  f(x,y) 


y(xn) 


we  compute  in  turn 


kl  =  h  f(Vyn> 

k2  =  h  f(xn  +  h/2,  y^  +  kx/2) 

k3  =  h  f(xn  +  h/2,  yn  +  k2/2) 

k.  =  h  f  (x  +  h, ,  y  +  k,) 

4       n    i*  Jn  3 

Vl  =  "n  +  i(kl  +  2k2  +  2k3  +  k4). 

This  process  may  be  described  in  geometric  terms.   At  the  point 

(x  ,  y  )  we  compute  the  slope  (k^  /h)  and  using  it  we  go  one  half  step 

forward  and  examine  the  slope  there.   Using  this  new  slope  (k»/h)  we 

again  start  at  (x  ,y  ),  go  one  half  step  forward,  and  again  sample  the 

slope.   Using  this  latest  slope  (k./h)  we  again  start  at  (x  ,y  )  but 

this  time  we  go  a  full  step  forward  where  we  examine  the  slope  (k,/h). 

The  four  slopes  are  averaged,  using  weights  1/6,  2/6,  2/6,  1/6,  and 

using  this  average  slope  we  make  the  final  step  from  (x  ,y  )  to  (x   ,, 
&  or  r      v  n*  n      n+1" 

y   ).   If  f(x,y)  did  not  depend  on  y,  then  the  averaging  would  produce 
Simpsons' s  formula.   The  method  has  an  error  term  proportional  to  h  , 

It  is  evident  that  the  method  throws  away  all  old  information  and 
begins  each  complete  step  anew,  and  hence  is  hardly  likely  to  be  as 
efficient  as  methods  which  take  advantage  of  old  information.   It  is 
also  evident  that  there  is  no  check  on  whether  the  step  size  is  too 
small  or  too  large,  though  perhaps  a  study  of  the  k.  might  give  a  clue; 
this  is  not  usually  done. 
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The  general  spirit  of  the  derivation  is  that  the  functions  f(x,y) 
which  are  on  the  right-hand  sides  are  all  expanded  in  series  in  powers  of 
h  and  the  corresponding  derivatives  are  equated  to  eliminate  the  lower 
powers  of  h. 

The  method  as  used  in  a  subroutine  for  the  investigations  of  this 
thesis  is  given  below  in  Fortran  language. 


SUBROUTINE  RKUTTA(N,T,X,DT) 

DIMENSION  X(300),  AK(4,300),  XD0T(300),  C(4) 

C(1)=0.0 

C(2)=0.5 

C(3)=0.5 

C(4)=1.0 

DO  4  1=1,4 

TC=T  +  C(I)*DT 

DO  2  J=1,N 

2  XC(J)=X(J)  +  C(I)*AK(I-1),J) 
CALL  DERIV  (TC,  XC,  XDOT,  N) 
DO  4  J=1,N 

4  AK(I,J)-DT*XDOT(J) 
DO   3  J=1,N 

3  X(J)+X(J)   +    (AK(1,J)   +2.*AK(3,J)    +  AK(4,J))/6. 
RETURN 

END 
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APPENDIX  Vb 
FORTRAN  LANGUAGE  PROGRAM  TO  OBTAIN  TRAJECTORIES  FOR  THE 
UNCONTROLLED  SYSTEM 
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•JOB  GIBSON    9  JAN  1963  SPECIAL  RUN 

PROGRAM  UNCON 

DIMENSION   X(30),   XH900),   YK900),   Y2(900),   Y3(900),   Y4(900) 
1,   Y5(900),   Y6(900),   Y7(900),   Y8(900),   Y9(900),   Y10(900), 
2Y1K900),   Y12(900) 

NUMPTS  =  0 

READ  3,  N,  TC,  TF,  DT,  (X(J),J=1,N) 

N  =  NO.  OF  EQNS.,   TO  ^INITIAL  VALUE.    TF  *  FINAL  VALUE 

DT  =  TIME  STEP.     X(J)  =  ARRAY  OF  DEPENDENT  VARIABLES 
3  FORMAT(I10/(8FlO.O)) 


104 


105 
20 


T  =  TO 

NUMPTS  = 


XI (NUMPTS) 
Yl(NUMPTS)  = 
Y2(NUMPTS)  = 
Y3(NUMPTS)  = 
Y4(NUMPTS)  = 
Y5(NUMPTS)  - 
Y6(NUMPTS)  = 
Y7(NUMPTS)  ■ 
Y8(NUMPTS)  = 
Y91NUMPTS)  = 
YIO(NUMPTS)  ■ 
Yll(NUMPTS)  = 
Y12(NUMPTS)  = 


NUMPTS  ♦  1 


T 

X(l  ) 

X(2) 

X(3) 

X(U) 

X(5) 

X(6) 

X(7) 

X{8) 

X(9) 
X{  10) 
X(ll) 
X(12) 


IF(.TF  -  T 
CALL  GRAPH 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


GRAPH 
GRAPH 
GRAPH 
GRAPH 
GRAPH 
GRAPH 
GRAPH 
GRAPH 
GRAPH 
GRAPH 
GRAPH 


DT)  10,20,20 
(NUMPTS, XI  ,Y12 
(NUMPTS, XI  ,Y11 
(NUMPTS, XI  ,Y1C 
(NUMPTS, XI,  Y9 
(NUMPTS, XI 
(NUMPTS, XI 
(NUMPTS, XI 
(NUMPTS, XI 
(NUMPTS, XI 
(NUMPTS, XI 
(NUMPTS, XI 


Y8 
Y7 
Y6 
Y5 
YU 
Y3 
Y2 


(NUMPTS,  *],  Yl 
FORMAT  (12H0TRAJECT0RY  / (6 
PRINT  104,  (  Yl  (I  ) 


8) 

8) 

8) 

8) 

8) 

8) 

8) 

8) 

8) 

8) 

8) 

8) 
16,9) ) 
,  NUMPTS) 
,  NUMPTS) 
,  NUMPTS) 
,  NUMPTS) 
,  NUMPTS) 
,  NUMPTS) 
,  NUMPTS) 
,  NUMPTS) 
,  NUMPTS) 
,  NUMPTS) 
,  NUMPTS) 
,  NUMPTS) 


PRINT  104,  (  Y2  (I) 

PRINT  104,  (  Y3  (I) 

PRINT  10U,  (  Y4  (I) 

PRINT  104,  (  Y5  (I) 

PRINT  104,  (  Y6  (I) 

PRINT  104,  {  Y7  (I) 

PRINT  104,  (  Y8  (I) 

PRINT  10U,  (  Y9  (I) 

PRINT  104,  (  Y10(I) 

PRINT  104,  (  Y1KI) 

PRINT  104,  (  Y12U) 

PRINT  105,  NUMPTS 

FORMAT  (10H0NUMPTS  ■    110) 

STOP 

CALL  RKUTTA(N,  T,  X,  DT) 

T  =  T  ♦  DT 

GO  TO  1 

END 

SUBROUTINE   RKUTTA( N, T,X ,DT ) 

DIMENSION  X(30),AK(4,30) , XDOT ( 30 ) , XC ( 30) ,C(4) 

C(l )=0.0 

C(2)=0.5 

C(3)=0.5 

C(4)=1.0 

DO  4  1=1,4 

TC=T+C( I )*DT 

DO  2  J  =  l  ,N 

XC(  J)  =X( J)      ♦  C( I)*AK( 1-1,  J) 

CALL  DERIV(TC,XC,XDOT) 

DO  4  J=1,N 


72 


AKU.J)  =  DT«XDCT(J) 

DO  3  J  =  I.N 

XU)=XCJ)*{AK(1  ,J)*2.»AK(2, J)+2.»AK(3,J)+AK(U, J)  )/6 

RETURN 

END 

SUBROUTINE  DERIV(T  ,X  ,X00T) 

DIMENSION  XDCT(30>, X(30) 

)+  SINE  (  X(  1) ) 
)+  SINF  {  X(  2) ) 
)+  SINF  (  Xt  3) ) 
)+  SINF  (  X{  4)) 
)+  SINF  (  X(  5) ) 
)+  SINF  (  X(  6)  ) 
)+  SINF  {  X(  7)) 


XDOT  (1) 

XDOT  (2) 

XDOT  (3) 

XDOT  U) 

XDOT  (5) 

XDOT  (6) 

XDOT  (7) 

XDOT  (8) 

XDOT  (9) 

XDOT  (  10  ) 

XDOT  { 

XDOT  (  12  ) 
END 
END 


(  1 
(  2 


0.1 

-0.1 
-0.1 
-0.1 
-C.l 
-0.1 
-0.1 
-C.l 

-0.1  »   X  (  9 
*  -0.1  »  X(  10  ) 
)  ~    -0.1  •  X(  11  )  + 
0.1  •  X{  12   ) 


(  8  )+  SINF  {  X(  8)  ) 
)+  SINF  (  X<  9)  ) 

SINF<  X(  10  )) 
SINF(  X(ll  })'' 
♦  SINFl  X(  12  )) 


12 


10.0 
•7.5 


0.1 
-8.0 


-.5 

-8.5 


X  VS  TIME  UNCONTROLLED  SYSTEM 
GIBSON  ANO  ALMSTEDT 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


-1.0       -2.0       -3. 
-10.0      -15.0      -20 
10    08 

7  DECEMBER  19.62 


-5.5 
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APPENDIX  Vc 

REPRESENTATIVE  FORTRAN  LANGUAGE  PROGRAM  TO  OBTAIN  COST 
FUNCTION  AS  A  FUNCTION  OF  INITIAL  P  (ADJOINT  VARIABLE) 
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.JOB  GIBSON    9  JAN  1963  SPECIAL  RUN   MAX  TIME  15  Mlff^^ 
PROG8AM  PJUX 

DIMENSION  X(950>,  XH950),  YK950) 
REAO  3,  TO,  TF,  DT,  XI,  AL,  AH,  DA 

TO=  INITIAL  VALUE,  TF=  FINAL  VALUE,  DT-  TIME  STEP,  XI*  INITIAL 
AH=  HIGHEST  INITIAL  ADJOINT,  AL=  LOWEST  INITIAL  AOJOINT, 
DA=  SPACING  CF  INITIAL  ADJCINTS  POINTS 
3  FORMAT  (8F10.C) 
N  =  900 

FORMS  MATRIX  OF  INITIAL  VALUES  OF  X 
DO  300  I  =  1 ,N,3 

300  X( I  )  =  XI 

FORMS  MATRIX  OF  P-ZERO 

X(2)  =  AL 

DO  301  I  *  5,N,3 

301  X(I )  =  X( 1-3)  ♦  DA 

ZEROIZES  INITIAL  VALUES  OF  COST  FUNCTION 
DO  302  I  *  3,  N,  3 

302  XU)  *  0.0 
T=  TO 

NUMPTS  =  N/3 
L  *  2 

DO  12  J  =  1, NUMPTS 
Xl(  J)  *  X(L) 

12  L  =  L«-3 

21  IF  (TF  -T>  10,  20,  20 

10  L  =  3 

DO  11  J  =  1, NUMPTS 
Yl( J)  =  X(L) 

11  L=  L+3 

DO  13   1=1,  NUMPTS 
J  =  2*1 1-1  )  ♦  I 

13  XU  )  =  X(  J) 
PRINT  101 

101  FORMAT  (  U3H0P-ZER0  VS    COST  FUNCTION  FINAL  X) 
PRINT  102  (  XKI),  YllI).  X(I),  1  =  1,  NUMPTS) 

102  FORMAT  (F8.2,  E16.9,   E16.9) 
CALL  GRAPH  (NUMPTS,  XI,  Yl,  8) 
STOP 

20  CALL  RKUTTA  (N,T,X,DT) 
T=  T+DT 
GO  TO  21 
END 

SUBROUTINE   RKUTTA( N,T ,X .DT ) 

DIMENSION  X(950),  AK(U,950>,  XD0T(950),  XC(950),  C(U) 
C( 1 )=0.0 
C(2)=0.5 
C(3)=0.5 
C(U)=1.0 


DO  k    1=1 .4 
TC=T+C(I)«DT 


DO  2  J»1,N 

2  XC(J)  =X( J)       +  C(I)»AKU-1.J) 

CALL  DERIV  ( TC ,  XC,  XDOT,  N) 
DO  U  J=1,N 
U  AK( I ,J)  «  DT«XDOT( J) 
DO  3  J  =  1  ,N 

3  X(J)=X(J)+(AK(1 ,J)+2.»AK(2, J)+2.»AK(3, J)+AK(U,J) J/6, 
RETURN 

END 

SUBROUTINE  DER I  V ( T ,X ,XDOT, N ) 
DIMENSION  XDCT(950),  X(950) 
DO  500  K=1,N,3 
L  =  K+1 
M  =  K  +  2 

XDOT(K}=  -0.1«  X(K)  +  SINF(X(K))+  .03125*  X(L) 
XDOT(L)=   0.1»  X(D-  X(D*  COSF(X(K))+  2.0»  X(K) 
500  XDOT(M)=  X(K>»»2  +  .015625«  X(L>*«2 
END 
END 


The  following  are  the  values  of  the  variables  read  in  for  this  program: 

TO  =  0.0 
TF  =  10.0 

DT  =  .1 
XI  =  20.0 
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AL  = 

-250.0 

AH  = 

50.0 

DA  = 

1.0 

APPENDIX  Vd 

REPRESENTATIVE  FORTRAN  LANGUAGE  PROGRAM  USED  TO  COMPUTE 
THE  OPTIMUM  CONTROL  POLICY  AND  TRAJECTORY  HAVING  THE 
P*(0)  FOR  THE  SYSTEM. 
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JOB  ALMSTEDT  MK2  MOO  1  1/20/63  MAX  TIME  5  MIN  X(0)  =  20.0  OPT.  TF(20) 
PROGRAM  MARK2  MOD  1 
7  DIMENSION  X(30),  X  1  (  10,201 ), U( 1 0,201 ) ,  CFCN(IO), 
71XRAY(201),  UNCLE1201).  TIME(201) 

Xl  =  TRAJECTORIES,  U=  CONTROL  EFFORT,  CFCN=  COST  FUNCTION 

READ  101,  N,  TO,  TF,  DT,  XI,  AL 

N=  TOTAL  NUMBER  OF  EONS,  TO  =  INITIAL  TIME.  TF  =  FINAL  TIME 

XI=  INITIAL  VALUE  OF  X,  AL  =  INITIAL  P  VALUE 

DO  1  I  =  1,N,3 

1  X( I )=  XI 

READ  102,  (X(I),  I»2fN,3J 

DO  5   1=  3,N,3 
5  X(I)=0.0 

T=  TO 

KK=  1 

NN  =  N/3 
30  IF(TF  -  T)  10,  20,  20 
To  L=  3 

DO  11  1=  1  ,  NN 

CFCN(I)=  X(L) 
11  L=L+3 

TIME( 1)  =  0.0 

DO  15  I  =  2,  KK 

15  TIME(I)=  TIME(I-1)+  DT 
DO  16   1=  1,  NN 
PRINT  103,  CFCN(I) 

DO  17  K  =  1 ,  KK 
XRAY  (K)  =  Xl( I,K) 
17  UNCLE  (K)  =  U(I,K) 

PRINT  104, (XRAY  (K),  K  =  1,  KK) 
PRINT  105,  (UNCLE  (K),  K=  1,KK) 
NPTS  =  KK 

CALL  GRAPH  (  NPTS,  TIME,  XRAY,  8) 
CALL  GRAPH  (  NPTS,  TIME,  UNCLE,  8) 

16  CALL  GRAPH  (  NPTS,  XRAY  ,  UNCLE,  8) 
STOP 

20  L=  1 

DO  21   1=  1,NN 
Xl( I  ,KK)  =  X(L) 
U(I,KK)  =  X(L+1)/8.0 

21  L  =  L+3 
KK=  KK*  1 

CALL  RKUTTA(N,T,  X,  DT) 
T=  T+  DT 
GO  TO  30 

101  FORMAT  (  110,  7F10.3) 

102  F0RMATC8F10.6) 

103  F0RMAT(17H0C0ST  FUNCTION  =  E16.9)) 

104  FORMAT( 12H0TRAJECT0RY  /(6E16.9) 

105  FORMAT( 16H0C0NTR0L  EFFORT  /(6E16.9)) 
END 

SUBROUTINE   RKUTTA (N, T,X , DT ) 

DIMENSION  X(601),  AK(4,601),  XD0T(601),  XC(601),  C(4) 

C(l )=0.0 

C(2)=0.5 

C(3)=0.5 

C(4)=1.0 

DO  U  1=1, 4 

TC=T+C(I)«DT 

DO  2  J=1,N 

2  XC(J)  =X(J)      +  C(I)*AK(I-1,J) 

CALL  DERIV  (TC,  XC,  XDOT,  N) 
DO  4  J=l ,N 
4  AMI  ,J)  =  DT»XDOT(J) 
DO  3  J  =  1,N 

3  X(  J)=X(  J)*(AK(1  ,J)+2.«AU2,J)+2.«AK(3,J)+AK(4,J)  )/6. 
RETURN 

END 

SUBROUTINE  DERI V ( T,X, XDOT, N) 

DIMENSION  XD0T(601),  X(601) 
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DO  500  K=1,N,3 
l=K+l 
M  =  K  +  2 

XDOT(K>*  -0,1*  X(K)*  SINF(X(K))+  .03125»  X(U 
XDOT(L)  *  O.U  XCLJ  -  X(L)«  COSFU(K))*  2.0*  X(K) 
500  XDOT(M)*  X(K)««2  ♦  .015625*  X(L)«»2 
ENO 
ENO 

90.0        20.0       0.1        20.0       0.0 
128.955976 

10    08  02 

2   TRAJECTORY   X  VS  TIME       P-2ER0  =  20.0  NO  1  TFC20) 
2   GIBSON  ANO   ALMSTEDT        JAN  1963   HK2  f 

o  i 

io   08  o; 

2   CONTROL  FUNCTION    U  VS  TIME        NO  1     TF(20) 
2   GIBSON  ANO   ALMSTEDT        JAN  1963   MK2 

0 

10    08  OJ 

2   INITIAL  CONTROL  VS.  TRAJECTORY   OPTIMAL  SYSTEM  TF(20) 
2   GIBSON  ANO  ALMSTEDT     JAN  1960   MK2 

0 
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APPENDIX  Ve 

FORTRAN  LANGUAGE  PROGRAM  UTILIZED  IN  THE  EVALUATION 
OF  THE  TRAJECTORIES  AND  CONTROL  POLICIES  FOR  THE 
SYSTEM  HAVING      J(X2)dt  AS  THE  COST  FUNCTION. 


79 


•.JOB  GIBSON        22  JANUARY  1963    MK4MCD0  ALT! 

PROGRAM   MK4  MCDO     ALT  1 
C      THIS  PROGRAM  USES  INTERGAL  X  SQUARE  AS  THE  COST  FUNCTION 

DIMENSION   X(30),  XR AY( 1 5 ,5C0 ) ,  U(15),  UNCLE ( 1 5, 500 ) ,  TIME(5Q0), 
1Y1 ( 101 ) ,  Y2( 101 ) ,  V( 15) 

READ  101, NN,  TO,  TF,  DT,  XI,  (U ( I ) , 1=1 ,NN ) 
101  FORMAT  (18/  15F5.0) 

N  =  2«  NN 
C      FORM   X-MATRIX 

DO  1  1=1, NN 

1  X(I)=XI 
N1=NN+1 

DO  2  1=  Nl ,N 

2  X(I )=  0.0 
T=  TO 
NPTS=  0 

27  IF(TF  -  T)  10,  20,  20 

20  NPTS  =  NPTS  +  1 

IF  (  1  -  NPTS)  21,26,26 

21  CALL  RKUTTA   (N,T,X,DT,  U,  V) 
T  =  T  +  DT 

26  TIME(NPTS)  =  T 
DO  28  I  =  1,  NN 
IF  (  X( I ) )  23,  24,  25 

23  UNCLE  (I ,NPTS)  =  -U( I ) 
GO  TO  22 

24  UNCLE  ( I, NPTS)  =  0.0 
GO  TO  22 

25  UNCLE  ( I,NPTS)=   U(I) 

22  XRAY  (I,NPTS)=  X(I) 

28  V(I)  =  UNCLE(I,NPTS) 
GO  TO  27 

10  DO  11  I  =  1, NN 
PRINT  111  ,  XU+NN) 

111  FORMAT  (17H0C0ST  FUNCTION  =  E16.9) 
PRINT  112,(XRAY  ( I , J ) , J= 1 , NPTS ) 

112  FORMAT  (12H0TRAJECT0RY  /  (6E16.9)) 
PRINT  113  (UNCLE  (I,J),J=  1,NPTS) 

113  FORMAT  (16H0CGNTROL  EFFORT  /(6E16.9)) 
DO  12  K=1,NPTS 

Yl  (K)  =  XRAY(I,K) 
12  Y2(K)  =  UNCLEf I,K) 

CALL  GRAPH  ( NPTS , TIME , Yl , 8 ) 

11  CALL  GRAPH  (NPTS, TIME, Y2, 8) 
PRINT  114,  NPTS 

114  FORMAT  {  8H0NPTS  =  15) 
STOP 

END 

SUBROUTINE    RKUTTA  (N,T,X,DT,  U,  V) 

DIMENSION   X(30),  AK(4,30),  XDOT  (30),  XC(30),  C(4),V(15),  UC30) 

C(1)=0.0 

C(2)=0.5 

C(3)=0.5 

C(4)=1.0 

DO  4  I  =  1 ,4 

TC  =  T  +  C(I )*DT 

DO  2  J  =1 ,N 

2  XC(J)  =  X(J)  +  C(I )*AK(I-1, J) 

CALL  DERIV  (  XC,  TC,-  XDOT,  N,  U,  V) 
DO  4  J  =  1,N 
4  AK( I ,J)  =  DT*XDCT( J) 
DO  3  J  =  1 ,N 

3  X( J)=X( J)+(AK( 1, J)+2.*AK(2, J)+2.*AK(3, J)+AK(4, J) )/6. 
RETURN 

END 

SUBROUTINE  DERIV  { X, T, XDOT, N,  U  ,  V) 

DIMENSION   XC0T(30),  X(30),  U(15),  V(15) 

NN  =  N/2 
Nl  =  NN  +  1 
DO  1.  1=  1  ,NN 

1  XDOT(I)  =  -.1*X(I)  +  SINF(Xd))  >.25*   V(I) 
DO  2  I  =  N1,N 

2  XDOT  (I)=   X(I-NN)»*2 
END 

END 

10 
0.0    10.0.1     20.0  -5.0-10.-15.  -20.  -25.  -30.  -35.  -40.  -45.  -50. 
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